Subroutine calsimeq(xx,pres,phi,totmol)

	USE ckvar

	IMPLICIT NONE
	REAL(dp) :: fj(nspc+1),matj(nspc+1,nspc+1),df(nspc+1)
	REAL(dp) :: xx(nspc+1),xold(nspc+1)
	REAL :: pres,phi
	REAL :: errmx,err(nspc+1),xxtot
	REAL :: totmol
	INTEGER :: i,j,k,n
	INTEGER :: nunk	! # of unknowns

	nunk=nspc+1
	xx(nspc+1)=totmol	! total # of moles is a unknown
!	do i=1,nspc
!	xx(i)=1.0/nspc
!	enddo

!	do 10 n=1,100
	xold=xx
	fj=0.0_dp
	matj=0.0_dp

! set up Kp equations for each reaction
! set up 1st derivation matrix of Kp equations
	do i=1,nreac
	write(*,*) 'Reaction#:',i,'Kp=',kp(i)
	fj(i)=1.0
	  do j=1,3
	  fj(i)=fj(i)*(xx(ispreac(i,j))*pres)**(nspreac(i,j))
	  enddo
	  do j=1,nunk
	    do k=1,3
	    if(j.eq.ispreac(i,k)) then
	    matj(i,j)=fj(i)*nspreac(i,k)/xx(j)
	    write(*,*) j,spc(ispreac(i,k)),xx(j),matj(i,i)
	    exit 
	    endif
	    enddo
!	  write(*,*) j,matj(i,j)
	  enddo
	fj(i)=fj(i)-kp(i)
	enddo

!! set up atom conservations
! 1. Sum of every mole fraction:Eq6
	do i=1,nspc
	fj(nreac+1)=fj(nreac+1)+xx(i)
	enddo
	fj(nreac+1)=fj(nreac+1)-1.0_dp
	matj(nreac+1,1)=1.0_dp
	matj(nreac+1,2)=1.0_dp
	matj(nreac+1,3)=1.0_dp
	matj(nreac+1,4)=1.0_dp
	matj(nreac+1,5)=1.0_dp
	matj(nreac+1,6)=1.0_dp
	matj(nreac+1,7)=1.0_dp
	matj(nreac+1,8)=1.0_dp
! 2. H elements conservation:Eq7
	fj(nreac+2)=2.0_dp*xx(1)+2.0_dp*xx(4)+xx(5)+xx(7)-2.0_dp/xx(9)
	matj(nreac+2,1)=2.0_dp
	matj(nreac+2,4)=2.0_dp
	matj(nreac+2,5)=1.0_dp
	matj(nreac+2,7)=1.0_dp
	matj(nreac+2,9)=2.0_dp/(xx(9)**2_dp)
! 3. O elements conservation:Eq8
	fj(nreac+3)=2.0_dp*xx(2)+xx(4)+xx(5)+xx(6)+xx(8)-1.0_dp/(phi*xx(9))
	matj(nreac+3,2)=2.0_dp
	matj(nreac+3,4)=1.0_dp
	matj(nreac+3,5)=1.0_dp
	matj(nreac+3,6)=1.0_dp
	matj(nreac+3,8)=1.0_dp
	matj(nreac+3,9)=1.0/(phi*(xx(9)**2.0_dp))
! 4. N elements conservation:Eq9
	fj(nreac+4)=2.0_dp*xx(3)+xx(8)-3.76_dp/(phi*xx(9))
	matj(nreac+4,3)=2.0_dp
	matj(nreac+4,8)=1.0_dp
	matj(nreac+4,9)=3.76_dp/(phi*(xx(9)**2.0_dp))

	CALL LU(matj,df,-fj)
!	write(*,*) 'df=',(df(i),i=1,9)	
	errmx=0.0
	do i=1,9
	xx(i)=xx(i)+df(i)
	err(i)=abs((xx(i)-xold(i))/xold(i))*100.0
	errmx=max(err(i),errmx)
	enddo
	write(*,*) n,'errmx=',errmx,(df(i),i=1,9)
	if(errmx.le.0.001) go to 20
!  10	continue
  20	xxtot=0.0_dp
	do i=1,9
	  xxtot=xxtot+xx(i)
	enddo
	totmol=xx(nspc+1)

End Subroutine

Subroutine LU(a,x,b)
! LU DECOMPOSITION MADE BY J.C. KIM :
	INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15,100)

	REAL(dp) :: a(9,9),x(9),b(9)
	REAL(dp) :: factor,sum
	integer :: n,i,j,k
	n=9

	do k=1,n-1
	  do i=k+1,n
	    factor=a(i,k)/a(k,k)
	    a(i,k)=factor
	    do j=k+1,n
	      a(i,j)=a(i,j) - factor*a(k,j)
	    enddo
	  enddo
	enddo

	do i=2,n
	  sum=b(i)
	  do j=1,i-1
	    sum=sum - a(i,j)*b(j)
	  enddo
	  b(i)=sum
	enddo

	x(n)=b(n)/a(n,n)
	do i=n-1,1,-1
	  sum=0.0
	  do j=i+1,n
	    sum=sum+a(i,j)*x(j)
	  enddo
	  x(i)=(b(i)-sum)/a(i,i)
	enddo

	return
end subroutine
